#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <time.h>

#define PI 3.141592653589793
#define G 4*pow(PI,2)
#define diretorioWindows "C:\\dados\\"

void jacobiana2C(double* j, double m1,double m2,double r1x,double r1y,double r2x,double r2y,int tamanho);
void abrirArquivos(char SO,char* metodo,double tFinal,double n);
void carregarArquivosLinha(double *z, int dimensao,int tamanho);
void F(double *J,double* delta,int tamanho,double *resultado);
double FY(double *deltaP,double* delta,double t,int tamanho);
double FW(double t,double y);
void produtoEscalarVetor(double* v,double a,int tamanho, double *resultado);
void somaVetores(double* v,double* v2,int tamanho, double *resultado);
void imprimir(double* v,int tamanho);
double calculaNorma(double *v,int dimensao);
void finalizar(void);

FILE *f[5];
double *z,*J,*R1,*R2,*R3,*R4,*R5,*R6,k1,k2,k3,k4,k5,k6,*delta,*deltaP,*temp,*temp2,*temp3,temp4,*norma,*y,*w;

int main(int argc, char *argv[])
{
	clock_t inicio=clock();

	double tFinal,n,dt;
	int i=0,j,dimensao,corpos,T,tamanho,normalizacao=10000;
	char aux[100];

	tFinal=atof(argv[1]);
	n=atof(argv[2]);

	dimensao=2;
	corpos=2;
	tamanho=2*corpos*dimensao;
	T=pow(2,n)+1;
	//T=2;
	dt=tFinal/pow(2,n);

	abrirArquivos('w',"estrelaDupla",tFinal,n);
	J=(double *)calloc(pow(tamanho,2), sizeof(double));
	delta=(double *)calloc(tamanho, sizeof(double));
	deltaP=(double *)calloc(tamanho, sizeof(double));
	temp=(double *)calloc(tamanho, sizeof(double));
	temp2=(double *)calloc(tamanho, sizeof(double));
	temp3=(double *)calloc(tamanho, sizeof(double));
	norma=(double*)calloc(T,sizeof(double));
	y=(double*)calloc(T,sizeof(double));
	w=(double*)calloc(T,sizeof(double));
	z=(double*)calloc(tamanho,sizeof(double));
	R1=(double*)calloc(tamanho,sizeof(double));
	R2=(double*)calloc(tamanho,sizeof(double));
	R3=(double*)calloc(tamanho,sizeof(double));
	R4=(double*)calloc(tamanho,sizeof(double));
	R5=(double*)calloc(tamanho,sizeof(double));
	R6=(double*)calloc(tamanho,sizeof(double));

	sprintf(aux,"%0.15lf ",w[0]);
	fputs(aux,f[4]);
	fputs("\n",f[4]);

	//delta inicial aleatorio
	srand( (unsigned)time(NULL) );
	for(i=0;i<tamanho;i++)
	{
		delta[i]=rand()/(double)RAND_MAX;
		delta[i]=delta[i]-0.5;
		//delta[i]=1;
	}

	// norma delta inicial nao normalizado
	norma[0]=calculaNorma(delta,tamanho);
	//normalizacao delta inicial
	for(i=0;i<tamanho;i++)
	{
		delta[i]=delta[i]/norma[0];
	}
	imprimir(delta,tamanho);

	//norma delta inicial normalizado
	norma[0]=calculaNorma(delta,tamanho);

	puts("");

	for(i=1;i<T;i++)
	{
		carregarArquivosLinha(z,dimensao,tamanho);
		jacobiana2C(J,1,1,z[0],z[1],z[2],z[3],tamanho);
		//calcula delta ponto no instante
		F(J,delta,tamanho,deltaP);

		//resolver y por RK45
		//c=[0, 1/4, 3/8, 12/13, 1, 1/2]';
		//b5=16/135, 0, 6656/12825, 28561/56430, -9/50, 2/55
		k1=FY(deltaP,delta,(i-1)*dt,tamanho)*dt;

		temp4=(i-1)*dt+dt*1.0/4;
		k2=FY(deltaP,delta,temp4,tamanho)*dt;

		temp4=(i-1)*dt+dt*3.0/8;
		k3=FY(deltaP,delta,temp4,tamanho)*dt;

		temp4=(i-1)*dt+dt*12.0/13;
		k4=FY(deltaP,delta,temp4,tamanho)*dt;

		temp4=(i-1)*dt+dt*1;
		k5=FY(deltaP,delta,temp4,tamanho)*dt;

		temp4=(i-1)*dt+dt/2;
		k6=FY(deltaP,delta,temp4,tamanho)*dt;

		temp4= 16.0/135*k1+6656.0/12825*k3+28561.0/56430*k4-9.0/50*k5+2.0/55*k6;
		y[i]=y[i-1]+temp4;

		//resolver w por RK45
		//c=[0, 1/4, 3/8, 12/13, 1, 1/2]';
		//b5=16/135, 0, 6656/12825, 28561/56430, -9/50, 2/55
		k1=FW(i*dt,y[i-1])*dt;

		temp4=(i-1)*dt+dt*1.0/4;
		k2=FW(temp4,y[i-1])*dt;

		temp4=(i-1)*dt+dt*3.0/8;
		k3=FW(temp4,y[i-1])*dt;

		temp4=(i-1)*dt+dt*12.0/13;
		k4=FW(temp4,y[i-1])*dt;

		temp4=(i-1)*dt+dt*1;
		k5=FW(temp4,y[i-1])*dt;

		temp4=(i-1)*dt+dt/2;
		k6=FW(temp4,y[i-1])*dt;

		temp4= 16.0/135*k1+6656.0/12825*k3+28561.0/56430*k4-9.0/50*k5+2.0/55*k6;
		w[i]=w[i-1]+temp4;

		//resolver delta por RK45
		F(J,delta,tamanho,temp);
		produtoEscalarVetor(temp,dt,tamanho,R1);
		produtoEscalarVetor(R1,0.25,tamanho,temp);
		somaVetores(temp,delta,tamanho,temp2);

		F(J,temp2,tamanho,temp);
		produtoEscalarVetor(temp,dt,tamanho,R2);
		produtoEscalarVetor(R1,3.0/32,tamanho,temp);
		somaVetores(temp,delta,tamanho,temp2);
		produtoEscalarVetor(R2,9.0/32,tamanho,temp);
		somaVetores(temp,temp2,tamanho,temp3);

		F(J,temp3,tamanho,temp);
		produtoEscalarVetor(temp,dt,tamanho,R3);
		produtoEscalarVetor(R1,1932.0/2197,tamanho,temp);
		somaVetores(temp,delta,tamanho,temp2);
		produtoEscalarVetor(R2,-7200.0/2197,tamanho,temp);
		somaVetores(temp,temp2,tamanho,temp3);
		produtoEscalarVetor(R3,7296.0/2197,tamanho,temp);
		somaVetores(temp,temp3,tamanho,temp2);

		F(J,temp2,tamanho,temp);
		produtoEscalarVetor(temp,dt,tamanho,R4);
		produtoEscalarVetor(R1,439.0/216,tamanho,temp);
		somaVetores(temp,delta,tamanho,temp2);
		produtoEscalarVetor(R2,-8.0,tamanho,temp);
		somaVetores(temp,temp2,tamanho,temp3);
		produtoEscalarVetor(R3,3680.0/513,tamanho,temp);
		somaVetores(temp,temp3,tamanho,temp2);
		produtoEscalarVetor(R4,-845.0/4104,tamanho,temp);
		somaVetores(temp,temp2,tamanho,temp3);

		F(J,temp3,tamanho,temp);
		produtoEscalarVetor(temp,dt,tamanho,R5);
		produtoEscalarVetor(R1,-8.0/27,tamanho,temp);
		somaVetores(temp,delta,tamanho,temp2);
		produtoEscalarVetor(R2,2.0,tamanho,temp);
		somaVetores(temp,temp2,tamanho,temp3);
		produtoEscalarVetor(R3,-3544.0/2565,tamanho,temp);
		somaVetores(temp,temp3,tamanho,temp2);
		produtoEscalarVetor(R4,1859.0/4104,tamanho,temp);
		somaVetores(temp,temp2,tamanho,temp3);
		produtoEscalarVetor(R5,-11.0/40,tamanho,temp);
		somaVetores(temp,temp3,tamanho,temp2);

		F(J,temp2,tamanho,temp);
		produtoEscalarVetor(temp,dt,tamanho,R6);
		produtoEscalarVetor(R1,16.0/135,tamanho,temp);
		produtoEscalarVetor(R2,0,tamanho,temp2);
		somaVetores(temp,temp2,tamanho,temp3);
		produtoEscalarVetor(R3,6656.0/12825,tamanho,temp);
		somaVetores(temp,temp3,tamanho,temp2);
		produtoEscalarVetor(R4,28561.0/56430,tamanho,temp);
		somaVetores(temp,temp2,tamanho,temp3);
		produtoEscalarVetor(R5,-9.0/50,tamanho,temp);
		somaVetores(temp,temp3,tamanho,temp2);
		produtoEscalarVetor(R6,2.0/55,tamanho,temp);
		somaVetores(temp,temp2,tamanho,temp3);
		somaVetores(delta,temp3,tamanho,delta);

		if(i%normalizacao==0)
		{
			norma[i]=calculaNorma(delta,tamanho);
			for(j=0;j<tamanho;j++)
			{
				delta[j]=delta[j]/norma[i];
			}
		}

		sprintf(aux,"%0.15lf ",w[i]/((i+1)*dt));
		fputs(aux,f[4]);
		fputs("\n",f[4]);
	}
	printf("%0.20lf\n",w[i-1]/((i)*dt));

	int minutos=0;
	int segundos=0;
	int horas=0;
	double tempo=0;
	int tempo2=0;

	clock_t fim=clock();
	clock_t total=fim-inicio;
	tempo=(double)total/CLOCKS_PER_SEC;
	tempo2=(int)tempo;
	horas=(tempo2/3600);
	minutos=(tempo2/60)-horas*60;
	segundos=(tempo2 % 60);
	printf("%d h %d min %d s\n",horas, minutos ,segundos);
	printf("%lf \n",tempo);

	void finalizar(void);
	return 0;
}

void imprimir(double* v,int tamanho)
{
	int i;
	for(i=0;i<tamanho;i++)
	{
		printf("%0.10lf\n",v[i]);
	}
	puts("");
}

void somaVetores(double* v,double* v2,int tamanho, double *resultado)
{
	int i;

	for(i=0;i<tamanho;i++)
	{
		resultado[i]=v[i]+v2[i];
	}
}

void produtoEscalarVetor(double* v,double a,int tamanho,double *resultado)
{
	int i;

	for(i=0;i<tamanho;i++)
	{
		resultado[i]=a*v[i];
	}
}

void F(double *J,double* delta,int tamanho, double *resultado)
{
	int i,j;
	double temp;

	for(i=0;i<tamanho;i++)
	{
		temp=J[i*tamanho]*delta[0];
		for(j=1;j<tamanho;j++)
		{
			temp+=J[i*tamanho+j]*delta[j];
		}
		resultado[i]=temp;
	}
}

double FY(double *deltaP,double* delta,double t,int tamanho)
{
	double temp1=0, temp2=0;
	int i;

	for(i=0;i<tamanho;i++)
	{
		temp1+=deltaP[i]*delta[i];
		temp2+=delta[i]*delta[i];
	}

	return (temp1/temp2)*t;
}

double FW(double t,double y)
{
	if(t==0)
	{
		t=1;
	}
	double temp=2*y/t;
	return temp;
}

void abrirArquivos(char SO,char* metodo,double tFinal,double n)
{
	char *barra,*diretorio,*temp,*raiz;

	raiz=(char*)calloc(1000,sizeof(char));
	diretorio=(char*)calloc(1000,sizeof(char));
	temp=(char*)calloc(1000,sizeof(char));

	if(SO=='w')
	{
		barra="\\";
	}
	else
	{
		barra="/";
	}

	strcpy(raiz,diretorioWindows);
	strcat(raiz,metodo);
	strcat(raiz,barra);
	sprintf(temp, "%0.0lf", tFinal);
	strcat(raiz,temp);
	strcat(raiz,barra);
	sprintf(temp, "%0.0lf", n);
	strcat(raiz,temp);
	strcat(raiz,barra);

	strcpy(diretorio,raiz);
	strcat(diretorio,"r1.txt");
	f[0]=fopen(diretorio,"r");

	strcpy(diretorio,raiz);
	strcat(diretorio,"r2.txt");
	f[1]=fopen(diretorio,"r");

	strcpy(diretorio,raiz);
	strcat(diretorio,"p1.txt");
	f[2]=fopen(diretorio,"r");

	strcpy(diretorio,raiz);
	strcat(diretorio,"p2.txt");
	f[3]=fopen(diretorio,"r");

	strcpy(diretorio,raiz);
	strcat(diretorio,"megno.txt");
	f[4]=fopen(diretorio,"w");

	free(raiz);
	free(diretorio);
	free(temp);
}

void carregarArquivosLinha(double *z, int dimensao,int tamanho)
{
	int k=0;
	char *temp,*temp2;

	for(k=0;k<tamanho/2;k++)
	{
		temp=(char*)calloc(1000,sizeof(char));
		if(!temp)
		{
			puts("carregar temp");
		}
		temp2=(char*)calloc(1000,sizeof(char));
		if(!temp2)
		{
			puts("carregar temp2");
		}
		fgets(temp,1000,f[k]);
		int i=0,j=0;
		do{
			temp2[i]=temp[i];
			i++;
		}while(temp[i]!=' ');
		z[2*k]=atof(temp2);
		do{
			temp2[j]=temp[i];
			i++;
			j++;
		}while(temp[i]!=' ');
		z[2*k+1]=atof(temp2);
		free(temp);
		free(temp2);
	}
}

void jacobiana2C(double* j, double m1,double m2,double r1x,double r1y,double r2x,double r2y,int tamanho)
{
	double G2=(double)G;
	double t100 = r1x-r2x;
	double t101 = r1y-r2y;
	double t106 = r1x*2;
	double t107 = r2x*2;
	double t102 = t106-t107;
	double t103 = pow(t100,2);
	double t104 = pow(t101,2);
	double t105 = t103+t104;
	double t108 = pow((1/t105),(5.0/2.0));
	double t109 = pow((1/t105),(3.0/2.0));
	double t110 = pow(t102,2.0);
	double t111 = G2 * m1 * m2 * t108 * t110* 3.0/4.0;
	double t112 = r1y*2;
	double t113 = r2y*2;
	double t114 = t112-t113;
	double t115 = G2*m1*m2*t102*t108*t114*3.0/4.0;
	double t116 = G2*m1*m2*t109;
	double t117 = pow(t114,2);
	double t118 = G2*m1*m2*t108*t117*(3.0/4.0);
	double t119 = -t111+t116;
	double t120 = t116-t118;
	double t121 = -t116+t118;

	int colunas=3*2*2;

	j[0*colunas+4]=1/m1;
	j[1*colunas+5]=1/m1;
	j[2*colunas+6]=1/m2;
	j[3*colunas+7]=1/m2;
	j[4*colunas+0]=t111-G2*m1*m2*t109;
	j[4*colunas+1]=t115;
	j[4*colunas+2]=t119;
	j[4*colunas+3]=-t115;
	j[5*colunas+0]=t115;
	j[5*colunas+1]=t121;
	j[5*colunas+2]=-t115;
	j[5*colunas+3]=t120;
	j[6*colunas+0]= t119;
	j[6*colunas+1]=-t115;
	j[6*colunas+2]=t111-t116;
	j[6*colunas+3]=t115;
	j[7*colunas+0]=-t115;
	j[7*colunas+1]=t120;
	j[7*colunas+2]=t115;
	j[7*colunas+3]=t121;

	//printf("%0.25lf\n",j[3*6+3]);

	/*int i,k;

	for(i=0;i<12;i++)
	{
		for(k=0;k<12;k++)
		{
			printf("%3.6lf ",j[i*colunas+k]);
		}
		puts("");
	}*/
}

double calculaNorma(double *v,int dimensao)
{
	double temp=0;
	int i;

	for(i=0;i<dimensao;i++)
	{
		temp+=v[i]*v[i];
	}
	temp=sqrt(temp);
	return temp;
}

void finalizar(void)
{
	fclose(f[0]);
	fclose(f[1]);
	fclose(f[2]);
	fclose(f[3]);
	fclose(f[4]);
	fclose(f[5]);
	fclose(f[6]);

	free(J);
	free(delta);
	free(temp);
	free(temp2);
	free(temp3);
	free(norma);
	free(y);
	free(w);
	free(z);
	free(R1);
	free(R2);
	free(R3);
	free(R4);
	free(R5);
	free(R6);
}


